Escaped Muons

1 Introduction

The goal of this document is to examine “escaped muons” that do not decay in the storage region and do not enter a calorimeter. Some of these muons may leave the storage ring altogether.

I generated art files of simulated events using mdc2 with different magnetic fields turned on or off (scenarios). We want to examine,

• Muons that are eaten by the iron yoke and muons that escape the world volume entirely.
• Where do muons escape the world (3D and heat map plots?)
• Comparison of the different scenarios

# Load necessary libraries
library(reticulate)  # Access to python
library(stringr)   # string functions
library(parallel)  # Parallel processing (built-in to R)
library(dplyr)     # data analysis
library(purrr)     # Functional programming
library(readr)     # Read formats
library(ggplot2)   # ploting
library(rgl)       # 3D plots
library(gridExtra) # Arranging plots
library(testthat)  # testing
#
library(readGallery) # Read art Gallery files

# Put all readGallery::useDataProduct calls here
useDataProduct('std::vector<gm2truth::TrackingActionArtRecord>')
useDataProduct('std::vector<gm2ringsim::GeantTrackRecord>')

2 Load the samples

Samples are located on the Fermilab dCache in the directory /pnfs/GM2/scratch/users/lyon/arr_20170321.

# Function to properly alter paths to accomodate XRootD
#    e.g. convert /pnfs/BLA --> root://fndca1.fnal.gov/pnfs/fnal.gov/usr/BLA
xrootd_ify <- function(aPath) paste0('root://fndca1.fnal.gov/pnfs/fnal.gov/usr', str_replace(aPath, '/pnfs', ''))

Determine the locations of the data files

# Determine locations of our files
system("ssh lyon@gm2gpvm04.fnal.gov 'ls /pnfs/GM2/scratch/users/lyon/arr_20170321/*/*.root'", intern=T) -> arrFiles
arrFiles %>% xrootd_ify() -> arrXrdFiles
arrXrdFiles

 [1] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root"   
 [2] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243526_0/ARR_unified_noDipole_cyl.root"     
 [3] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243527_0/ARR_unified_noInflector_cyl.root"  
 [4] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243528_0/ARR_unified_noKickerQuads_cyl.root"
 [5] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243530_0/ARR_unified_noKicker_cyl.root"     
 [6] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243532_0/ARR_unified_noQuads_cyl.root"      
 [7] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243540_0/ARR_unified_everything_cyl.root"   
 [8] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243542_0/ARR_unified_noDipole_cyl.root"     
 [9] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243544_0/ARR_unified_noInflector_cyl.root"  
[10] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243547_0/ARR_unified_noKickerQuads_cyl.root"
[11] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243548_0/ARR_unified_noKicker_cyl.root"     
[12] "root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243552_0/ARR_unified_noQuads_cyl.root"      

Pull out the scenaios (these will be in the same order as the file names)

scenarios <- str_match(arrFiles, 'arr_20170321/.+unified_(.+)_cyl')[,2]
scenarios %>% unique()

[1] "everything"    "noDipole"      "noInflector"   "noKickerQuads" "noKicker"      "noQuads"      

Let’s look at the contents of these files (they’re all the same, so we’ll just do one) gm2 v7_04_00 will have product_sizes_dumper.

stringToRun <- str_interp(
  "ssh lyon@gm2gpvm04.fnal.gov 'source /cvmfs/gm2.opensciencegrid.org/prod7/g-2/setup
  setup gm2 v7_04_00 -q prof ; 
  product_sizes_dumper ${aFile} | grep ::'", 
  list(aFile = arrFiles[1])
)
#
system(stringToRun)

        1410762186      0.456 gm2truth::MagnetIronAndCryoArtRecords_artg4_magnetIronAndCryostat_mdc2arr.
         660570363      0.214 gm2truth::TrackingActionArtRecords_artg4__mdc2arr.
         637307490      0.206 gm2truth::TrajectoryArtRecords_artg4__mdc2arr.
         240738334      0.078 gm2ringsim::GeantTrackRecords_artg4_KeepStepsAction_mdc2arr.
         118035827      0.038 gm2truth::RingArtRecords_artg4_Ring_mdc2arr.
          12078805      0.004 gm2truth::GhostDetectorArtRecords_artg4_GhostCylinderDetector_mdc2arr.
           4482185      0.001 art::RNGsnapshots_randomsaver__mdc2arr.
           4086898      0.001 gm2truth::RingTrackingPlaneArtRecords_artg4_RingTrackingPlanes_mdc2arr.
           1373178      0.000 gm2truth::GhostDetectorArtRecords_artg4_GhostNearWorldDetector_mdc2arr.
           1185772      0.000 gm2truth::PhaseSpaceArtRecord_artg4__mdc2arr.
             36971      0.000 art::TriggerResults_TriggerResults__mdc2arr.

3 Escaping muons from the Tracking Action data

The tracking action data has birth and death data for every Geant track.

3.1 Load Tracking Action data

We’ll capture this information for the primary muon and its first generation daughters. See gm2dataproducts/mc/actions/track/TrackingActionArtRecord.hh.

We need a readGallery reader python class. Generate with,

readerClassSkel('gm2truth::TrackingActionArtRecord', writeFile = 'trackingActionReader.py')

Here is the reader,

read_file('trackingActionReader.py') %>% cat

from readGallery import GalleryReaderBase  # Necessary for the base class

class TrackingActionArtRecordReader(GalleryReaderBase):
  def __init__(self, inputTag):
    GalleryReaderBase.__init__(self, inputTag)
    self.names = ['fileEntry', 'eventEntry', 'trackType', 'trackID', 
                  'parentTrackID', 'turn', 'volumeUID', 'status',  'p', 'e', 'x', 'y', 'z', 'px', 'py', 'pz']

  def prepare(self, ROOT, ev):
    GalleryReaderBase.prepare(self, ROOT, ev)
    self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2truth.TrackingActionArtRecord))

  def fill(self, ROOT, ev):

    validHandle = self.getValidHandle(self.inputTag)  # Get the valid handle for 
                                                      # gm2truth::TrackingActionArtRecord

    if not validHandle.empty():                       # Does it have data?

      p = validHandle.product()                       # Get the corresponding data product (maybe a vector)

      # Fill from gm2truth::TrackingActionArtRecord
      for e in p:                             # Loop over elements and fill
        
        # Only accept the primary muon death and its first generation daughters birth
        if (e.trackID == 1 and e.status == 1) or (e.parentTrackID == 1 and e.status == 0):
        
          self.vals.append(
            [ ev.fileEntry(), ev.eventEntry(), e.trackType, e.trackID, 
              e.parentTrackID, e.turn, e.volumeUID, e.status, e.p, e.e, 
              e.x, e.y, e.z, e.px, e.py, e.pz ])

    return True

Make the reader class and object

trackingActionReaderC <- createReaderClass_from_file('trackingActionReader.py')$TrackingActionArtRecordReader

We have many files to process. Let’s do this reading in parallel.

We need to create a reader object per parallel job

trackingActionReaders <- map(1:length(arrXrdFiles), function(i) trackingActionReaderC(artInputTag("artg4")))

Load the data

jobs <- lapply(1:length(arrXrdFiles), function(i) getGalleryData(arrXrdFiles[i], trackingActionReaders[[i]]) %>%
                                                    mcparallel(name=i))

Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root
Closing file, read 663158402 bytes in 161 transactions
Timings: Overall time = 200.338 s
Time per event: min=0.000   mean=0.019   max=7.364
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243526_0/ARR_unified_noDipole_cyl.root
Closing file, read 657073380 bytes in 158 transactions
Timings: Overall time = 228.757 s
Time per event: min=0.000   mean=0.022   max=27.028
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243527_0/ARR_unified_noInflector_cyl.root
Closing file, read 668231087 bytes in 160 transactions
Timings: Overall time = 208.320 s
Time per event: min=0.000   mean=0.020   max=16.990
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243528_0/ARR_unified_noKickerQuads_cyl.root
Closing file, read 664723747 bytes in 161 transactions
Timings: Overall time = 200.955 s
Time per event: min=0.000   mean=0.019   max=8.426
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243530_0/ARR_unified_noKicker_cyl.root
Closing file, read 666316280 bytes in 160 transactions
Timings: Overall time = 234.668 s
Time per event: min=0.000   mean=0.022   max=17.103
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243532_0/ARR_unified_noQuads_cyl.root
Closing file, read 664218464 bytes in 136 transactions
Timings: Overall time = 535.363 s
Time per event: min=0.000   mean=0.053   max=334.559
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243540_0/ARR_unified_everything_cyl.root
Closing file, read 660527554 bytes in 161 transactions
Timings: Overall time = 266.787 s
Time per event: min=0.000   mean=0.026   max=21.746
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243542_0/ARR_unified_noDipole_cyl.root
Closing file, read 656304716 bytes in 158 transactions
Timings: Overall time = 209.096 s
Time per event: min=0.000   mean=0.021   max=17.045
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243544_0/ARR_unified_noInflector_cyl.root
Closing file, read 668176375 bytes in 160 transactions
Timings: Overall time = 209.131 s
Time per event: min=0.000   mean=0.020   max=11.414
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243547_0/ARR_unified_noKickerQuads_cyl.root
Closing file, read 664807569 bytes in 161 transactions
Timings: Overall time = 216.531 s
Time per event: min=0.000   mean=0.021   max=17.764
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243548_0/ARR_unified_noKicker_cyl.root
Closing file, read 666316217 bytes in 160 transactions
Timings: Overall time = 227.132 s
Time per event: min=0.000   mean=0.022   max=13.152
Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243552_0/ARR_unified_noQuads_cyl.root
Closing file, read 665617262 bytes in 161 transactions
Timings: Overall time = 221.986 s
Time per event: min=0.000   mean=0.021   max=9.361

trackingActionReturns <- mccollect(jobs)

Let’s merge all of the output.

trackingActionDF <- map_df(1:length(trackingActionReaders), 
                           function(i) galleryReader_df(trackingActionReaders[[i]]) %>% mutate(scenario=scenarios[i], 
                                                                                               fileEntry=i))
trackingActionDF

write_rds(trackingActionDF, 'trackingActionDF.rds')

How many muons per file?

trackingActionDF %>% filter(trackID==1) %>% group_by(scenario, fileEntry) %>% tally()

3.1.1 Add Volume Names

We need to turn the volume IDs into volume names. The volume IDs change for each file. We can run a small art job to determine the volume ID to name tables.

runForVolIDString <- function(i) {
  str_interp(
    "PVS_CSVOUT=${csvout}_${i}_volNames.csv gm2 -c ${fclPath}/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 ${aFile}",
    list( csvout=scenarios[i], i=i, fclPath=Sys.getenv("MRB_BUILDDIR"), aFile=arrXrdFiles[i]) 
  )
}
runForVolIDStrings <- sapply(1:length(arrXrdFiles), runForVolIDString)
runForVolIDStrings

 [1] "PVS_CSVOUT=everything_1_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root"       
 [2] "PVS_CSVOUT=noDipole_2_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243526_0/ARR_unified_noDipole_cyl.root"           
 [3] "PVS_CSVOUT=noInflector_3_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243527_0/ARR_unified_noInflector_cyl.root"     
 [4] "PVS_CSVOUT=noKickerQuads_4_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243528_0/ARR_unified_noKickerQuads_cyl.root" 
 [5] "PVS_CSVOUT=noKicker_5_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243530_0/ARR_unified_noKicker_cyl.root"           
 [6] "PVS_CSVOUT=noQuads_6_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243532_0/ARR_unified_noQuads_cyl.root"             
 [7] "PVS_CSVOUT=everything_7_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243540_0/ARR_unified_everything_cyl.root"       
 [8] "PVS_CSVOUT=noDipole_8_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243542_0/ARR_unified_noDipole_cyl.root"           
 [9] "PVS_CSVOUT=noInflector_9_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243544_0/ARR_unified_noInflector_cyl.root"     
[10] "PVS_CSVOUT=noKickerQuads_10_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243547_0/ARR_unified_noKickerQuads_cyl.root"
[11] "PVS_CSVOUT=noKicker_11_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243548_0/ARR_unified_noKicker_cyl.root"          
[12] "PVS_CSVOUT=noQuads_12_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243552_0/ARR_unified_noQuads_cyl.root"            

Let’s run the art jobs in parallel

jobs <- lapply(1:length(arrXrdFiles), function(i) mcparallel(system( runForVolIDStrings[i], intern=T ), name=i))
res <- mccollect(jobs)

running command 'PVS_CSVOUT=noKickerQuads_10_volNames.csv gm2 -c /home/me/gm2/renee_arr/build_slf7.x86_64/gm2analyses/fcl/physicalVolumeStoreToFile.fcl -n 1 root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243547_0/ARR_unified_noKickerQuads_cyl.root' had status 20

and load the CSV files.

jobs <- lapply(1:length(arrFiles), function(i) 
  str_interp("${scn}_${i}_volNames.csv", list(scn=scenarios[i], i=i)) %>%
    read_csv(col_names=c("volumeUID", "volName")) %>% mcparallel(name=i) )
#
volNameTables <- mccollect(jobs)
volNameDF <- map_df(1:length(volNameTables), function(i) 
                      volNameTables[[i]] %>% mutate(scenario=scenarios[i],
                                                    fileEntry=i))
volNameDF

write_rds(volNameDF, 'volNameDF.rds')

Now we merge with the tracking action data frame

trackingActionDF %>% inner_join(volNameDF) %>% 
  select(scenario, fileEntry, eventEntry, volName, everything()) -> trackingActionNDF

Joining, by = c("fileEntry", "volumeUID", "scenario")

trackingActionNDF

3.2 How many muons decay in the storage region?

These muons decay in the storage region (not quite the same as “stored”),

trackingActionNDF %>% filter(trackID==1, volName=='ArcSection[00]') %>% 
  group_by(scenario, fileEntry) %>% tally()

Interesting! So, from this, only 0.575% decay in the storage region for the “all on” (or everything) sample. That seems rather low. But continuing.

write_rds(trackingActionNDF, 'trackingActionNDF.rds')

3.3 Where do muons go to die?

trackingActionNDF %>% filter(trackID == 1) %>% group_by(scenario, volName) %>% tally()

It may be nice to categorize the volume names by area (e.g. “Arc”, “bellows”, “Ring”, “Inflector”).

# Pull out the first word in camel case (see http://stackoverflow.com/questions/29916065/how-to-do-camelcase-split-in-python)
camelCaseSplit <- function(s) str_split(s, '(?<=[a-z])(?=[A-Z])|(?<=[A-Z])(?=[A-Z][a-z])')

# Test the above
expect_equal(camelCaseSplit("RingPoleTipUpper")[[1]], c("Ring", "Pole", "Tip", "Upper"))
expect_equal(camelCaseSplit("ringPoleTipUpper")[[1]], c("ring", "Pole", "Tip", "Upper"))
expect_equal(camelCaseSplit("world")[[1]], c("world"))
expect_equal(camelCaseSplit("ringPole")[[1]], c("ring", "Pole"))

# We really just want the first word
takeFirstCamel <- function(s) camelCaseSplit(s) %>% map_chr(function(s) s[1])
takeFirstCamel("RingPoleTipUpper") %>% expect_equal("Ring")
takeFirstCamel("ringPoleTipUpper") %>% expect_equal("ring")
takeFirstCamel("world") %>% expect_equal("world")
takeFirstCamel("ringPole") %>% expect_equal("ring")
takeFirstCamel(c("ringPoleTipUpper", "world")) %>% expect_equal(c("ring", "world"))

Don’t confuse RingTrackingPlane with the iron

makeVolCategory <- function(s) ifelse( str_detect(s, 'RingTracking'),
                                       "RngTrkPlane",
                                       takeFirstCamel(s))

makeVolCategory(c("RingBottom", "RingTrackingPlane[8]", "world")) %>% expect_equal(c("Ring", "RngTrkPlane", "world"))

trackingActionNDF %>% mutate(volCategory = volName %>% makeVolCategory()) -> trackingActionNVDF
trackingActionNVDF

Show count by category

trackingActionNVDF %>% filter(trackID == 1) %>% group_by(scenario, volCategory) %>% tally() -> volCatTally
volCatTally %>% mutate(perc = n/20000 * 100) -> volCatTally
volCatTally

Let’s plot!

trackingActionNVDF %>% filter(trackID ==1) %>% ggplot(aes(x=volCategory, group=scenario)) + 
  geom_bar()

ggplot(volCatTally, aes(x = volCategory, y=perc, fill=scenario)) + 
  geom_bar(stat="Identity", width=0.5, position="dodge")

Can’t really see non-Ring/world ones. Let’s split them out.

# Function to make this easy
plotit <- function(d) ggplot(d, aes(x = volCategory, y=perc, fill=scenario) ) + 
  geom_bar(stat="Identity", width=0.5, position="dodge") + 
  xlab("Volume category") + ylab("Percent")
volCatTally %>% filter(volCategory %in% c("Ring", "world")) %>% plotit -> p1
volCatTally %>% filter(! volCategory %in% c("Ring", "world")) %>% plotit -> p2
grid.arrange(p1, p2)

Let’s see where they die

# See https://rpubs.com/hadley/97970 for how to wrap a multipart ggplot2 component
plotCommon <- function () {
  list(
    facet_wrap(~scenario),
    guides(col = guide_legend(override.aes = list(size=5, alpha=1))),
    scale_color_discrete(name="Volume\nCategory"),
    labs(x="z (mm)", y="x (mm)"),
    theme_minimal()
  )
}

trackingActionNVDF %>% filter(trackID == 1) %>% 
  ggplot( aes(x=z, y=x, color=volCategory)) + 
    geom_point(alpha=0.1) + 
    ggtitle('Where Muons Die') +
    plotCommon()

Here it is without the ring losses, which dominate

trackingActionNVDF %>% filter(trackID == 1, volCategory != 'Ring') %>% 
  ggplot( aes(x=z, y=x, color=volCategory)) + 
    geom_point(alpha=0.5) + 
    ggtitle("Where Muons Die", subtitle = "Ring escapees excluded for clarity") +
    plotCommon()

3.4 Does it make sense for so many muons to stop in the iron?

We have 3 GeV muons. Does it make sense that the vast majority of them stop in the iron? A table shows the Continuous Slow Down Approximation range of muons (CSDA). For a 3 GeV muon in iron, the CSDA is \(1.825 \times 10^{3}\) g/cm\(^2\). The density \(\rho\) is 7.874 g/cm\(^3\). So, where \(R\) is range,

\[R = \text{CSDA}/\rho\] Range is thus 2.32 m. This is begger than the width of the iron, but remember that the muons come in at an oblique angle. We’ll have to prove this.

4 How far do muons go in the iron?

We have the Geant stepping information, so we should be able to figure out how far the muons go in the iron.

Here is the Gallery reader,

read_file('GeantTrackReader.py') %>% cat

from readGallery import GalleryReaderBase  # Necessary for the base class

class GeantStepReader(GalleryReaderBase):
  def __init__(self, inputTag, eventsToFill = []):
    GalleryReaderBase.__init__(self, inputTag)
    self.names = ['fileEntry', 'eventEntry', 'trackID', 'globalStepNum', 'totalEnergyDeposit', 
                  'stepLength', 'volumeUID', 
                  'x', 'y', 'z', 'px', 'py', 'pz']
    self.eventsToFill = eventsToFill
    print self.eventsToFill

  def prepare(self, ROOT, ev):
    GalleryReaderBase.prepare(self, ROOT, ev)
    self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2ringsim.GeantTrackRecord))

  def fill(self, ROOT, ev):

    # Do we care about this event - if not then don't bother looking?
    if self.eventsToFill and not ev.eventEntry() in self.eventsToFill:
      return True
      
    print 'Reading entry %s' % ev.eventEntry()

    validHandle = self.getValidHandle(self.inputTag)  # Get the valid handle for gm2ringsim::GeantTrackRecord

    if not validHandle.empty():                       # Does it have data?

      p = validHandle.product()                       # Get the corresponding data product (maybe a vector)

      # Fill from gm2ringsim::GeantStep
      for e in p:                             # Loop over elements and fill
      
        # Get the steps
        for f in e.steps():

          self.vals.append(
            [ ev.fileEntry(), ev.eventEntry(), e.trackID(), f.globalStepNum(), 
            f.totalEnergyDeposit(), f.stepLength(), f.volumeUID(), 
            f.pos()[0], f.pos()[1], f.pos()[2], 
            f.p()[0], f.p()[1], f.p()[2] ])
            
        # If we're on the last event, don't bother reading any more
        if self.eventsToFill and ev.eventEntry() == self.eventsToFill[-1]:
          return False

    return True

geantStepReaderC <- createReaderClass_from_file('GeantTrackReader.py')$GeantStepReader

4.1 Muon distance in iron for one event

Let’s find an event where the muon did in the iron.

trackingActionNVDF %>% filter(trackID == 1, volName == 'RingYokeBottom')

Event 4 looks good.

geantStepOne <- geantStepReaderC(artInputTag('artg4:KeepStepsAction'), tuple(as.integer(4)))

(4,)

getGalleryData(arrXrdFiles[1], geantStepOne)

Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root
Reading entry 4

gsOneDF <- galleryReader_df(geantStepOne)

Merge in volume names

gsOneDF %>% mutate(fileEntry = 1) %>%
  inner_join(volNameDF) %>% 
  select(volName, everything()) -> gsOneDF

Joining, by = c("fileEntry", "volumeUID")

gsOneDF

Let’s make a 3D plot of this path

ring <- cylinder3d( center=rbind(c(0,-90, 0), c(0,90,0)),
                    radius=7112,
                    sides=20, closed = F)

Let’s do volume categories again.

gsOneDF %>% mutate(volCategory = volName %>% makeVolCategory %>% as.factor) -> gsOneDF

Let’s plot the path. Note that I have to make the legend separately since legend3d looks awful.

# lengend3d looks terrible - so do a regular legend
plot(1, type='n', axes=FALSE, ann=FALSE)
with(gsOneDF,
  legend("top", legend=unique(volCategory), col=as.integer(unique(volCategory)), pch=19)
)

clear3d()
with(gsOneDF,{
     plot3d(x=x, y=y, z=z, type='p', col = as.integer(volCategory), ylim=c(-400, 400))
})
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()

How much energy loss and distance traveled in each volume category

gsOneDF %>% group_by(volCategory) %>% summarize(totalDist_meters=round(sum(stepLength)/1000, 2), 
                                                totalELoss_MeV=round(sum(totalEnergyDeposit), 1),
                                                nSteps = n())

Ok - I believe this (note that the y scale is very small compared to the x,z scales - so the muon acutally travels quite far in x & z). The muon travels about 2.4m in the iron and ranges out. Physics (and Geant) work!

Let’s try another one

geantStepOneA <- geantStepReaderC(artInputTag('artg4:KeepStepsAction'), tuple(as.integer(10)))

(10,)

getGalleryData(arrXrdFiles[1], geantStepOneA)

Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root
Reading entry 10
Timings: Overall time = 178.569 s
Time per event: min=0.000   mean=0.100   max=1.101

gsOneADF <- galleryReader_df(geantStepOneA)
gsOneADF %>% mutate(fileEntry = 1) %>%
  inner_join(volNameDF) %>% 
  mutate(volCategory = volName %>% makeVolCategory %>% as.factor)%>% 
  select(volName, volCategory, everything()) -> gsOneADF

Joining, by = c("fileEntry", "volumeUID")

gsOneADF

gsOneADF %>% group_by(volCategory) %>% summarize(totalDist_meters=round(sum(stepLength)/1000, 2), 
                                                totalELoss_MeV=round(sum(totalEnergyDeposit), 1),
                                                nSteps = n()) -> gsOneADFSum
gsOneADFSum

Check sums

gsOneADFSum %>% summarize(sum(totalELoss_MeV))

gsOneADF %>% summarize(sum(totalEnergyDeposit))

# lengend3d looks terrible - so do a regular legend
plot(1, type='n', axes=FALSE, ann=FALSE)
with(gsOneADF,
  legend("top", legend=unique(volCategory), col=as.integer(unique(volCategory)), pch=19)
)

clear3d()
with(gsOneADF,{
     plot3d(x=x, y=y, z=z, type='p', col = as.integer(volCategory), ylim=c(-400, 400), xlim=c(-10000, 10000),
            zlim=c(-10000, 10000))
})
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()

This one scatters in the inflector.

4.2 Muon distance for many events

Can we show how far muons go in various materials?

Let’s collect data for 100 events.

We have a reader that can do n events.

read_file('GeantTrackReaderN.py') %>% cat

from readGallery import GalleryReaderBase  # Necessary for the base class

class GeantStepReaderN(GalleryReaderBase):
  def __init__(self, inputTag, startEvent=100, nEvents=10):
    GalleryReaderBase.__init__(self, inputTag)
    self.names = ['fileEntry', 'eventEntry', 'trackID', 'globalStepNum', 'globalTime', 'totalEnergyDeposit', 
                  'stepLength', 'volumeUID', 
                  'x', 'y', 'z', 'px', 'py', 'pz']
    self.startEvent = startEvent
    self.endEvent   = startEvent + nEvents
    
  def prepare(self, ROOT, ev):
    GalleryReaderBase.prepare(self, ROOT, ev)
    self.getValidHandle = ev.getValidHandle(ROOT.vector(ROOT.gm2ringsim.GeantTrackRecord))

  def fill(self, ROOT, ev):

    # Do we care about this event - if not then don't bother looking
    if ev.eventEntry() < self.startEvent:
      return True
      
    if ev.eventEntry() >= self.endEvent:
      return False

    validHandle = self.getValidHandle(self.inputTag)  # Get the valid handle for gm2ringsim::GeantTrackRecord

    if not validHandle.empty():                       # Does it have data?

      p = validHandle.product()                       # Get the corresponding data product (maybe a vector)

      # Fill from gm2ringsim::GeantStep
      for e in p:                             # Loop over elements and fill
      
        # Get the steps
        for f in e.steps():

          self.vals.append(
            [ ev.fileEntry(), ev.eventEntry(), e.trackID(), f.globalStepNum(), f.globalTime(), 
            f.totalEnergyDeposit(), f.stepLength(), f.volumeUID(), 
            f.pos()[0], f.pos()[1], f.pos()[2], 
            f.p()[0], f.p()[1], f.p()[2] ])

    return True

geantStepReaderNC <- createReaderClass_from_file('GeantTrackReaderN.py')$GeantStepReaderN

geantStepN <- geantStepReaderNC(artInputTag('artg4:KeepStepsAction'), 1000, 1000)

getGalleryData(arrXrdFiles[1], geantStepN)

Opening first file...
Successfully opened file root://fndca1.fnal.gov/pnfs/fnal.gov/usr/GM2/scratch/users/lyon/arr_20170321/15243525_0/ARR_unified_everything_cyl.root
Timings: Overall time = 66.083 s
Time per event: min=0.000   mean=0.030   max=2.367

gsNDF <- galleryReader_df(geantStepN)
gsNDF %>% mutate(fileEntry = 1,
                 p = sqrt(px*px + py*py + pz*pz)) %>%
  inner_join(volNameDF) %>% 
  mutate(volCategory = volName %>% makeVolCategory %>% as.factor)%>% 
  select(volName, volCategory, p, everything()) -> gsNDF

Joining, by = c("fileEntry", "volumeUID")

gsNDF

Let’s do a distribution of distance in the different volumes

gsNDF %>% group_by(eventEntry, volCategory) %>% summarize(
                                                totalDist_meters=round(sum(stepLength)/1000, 2), 
                                                totalELoss_MeV=round(sum(totalEnergyDeposit), 1),
                                                nSteps = n()) -> gsNDFSum
gsNDFSum

gsNDFSum %>% filter(volCategory == 'Ring') %>% 
  ggplot(aes(x=totalDist_meters)) + geom_histogram(bins=50) + 
    labs(x = 'Total distance in Ring material (m)', y='count', title='Muon in Ring Material')

Let’s look at this across the volume categories for energy loss

gsNDFSum %>% ggplot(aes(x=totalELoss_MeV)) + geom_histogram(bins=50) +
  facet_wrap(~volCategory, scales = "free") + 
  labs(x = 'Total energy loss (MeV)', y='count', title='Energy loss in material',
       subtitle = 'Note the different scales')

Let’s get an energy loss curve

gsNDF %>% ggplot(aes(x=p, y=totalEnergyDeposit)) + geom_point() + facet_wrap(~volCategory) + 
  labs(x="Muon Momentum (MeV/c)", y="Total Energy loss (MeV))")

Well, not sure what all this means. Let’s look for a particular event…

gsNDF %>% filter(eventEntry == 1300) %>% ggplot(aes(x=p, y=totalEnergyDeposit)) + geom_point() + facet_wrap(~volCategory) + 
  labs(x="Muon Momentum (MeV/c)", y="Total Energy loss (MeV))", title='Energy loss vs. Momentum for muons in event 1300' )

How many steps do we take?

gsNDFSum %>% ggplot(aes(x=nSteps)) + geom_histogram(bins=50) +
  facet_wrap(~volCategory, scales = "free") + 
  labs(x = 'Number of steps', y='count', title='Number of steps in materials',
       subtitle = 'Note the different scales')

Let’s look at two events

4.2.1 Muon with long distance in iron

gsNDFSum %>% filter(volCategory == 'Ring', totalDist_meters > 3)

Let’s look at event 1758

gsNDF %>% filter(eventEntry == 1758) -> gsNDF1758
gsNDF1758

gsNDF1758 %>% group_by(eventEntry, volCategory) %>% summarize(
                                                totalDist_meters=round(sum(stepLength)/1000, 2), 
                                                totalELoss_MeV=round(sum(totalEnergyDeposit), 1),
                                                nSteps = n()) -> gsNDF1758Sum
gsNDF1758Sum

Let’s plot it!

gsNDF1758 %>% mutate(volCategory = as.factor(as.character(volCategory))) -> gsNDF1758

# lengend3d looks terrible - so do a regular legend
plot(0, type='n', axes=FALSE, ann=FALSE)
with(gsNDF1758,
  legend("bottom", legend=unique(volCategory), col=as.integer(unique(volCategory)), pch=19)
)

clear3d()
with(gsNDF1758,{
     plot3d(x=x, y=y, z=z, type='p', col = as.integer(volCategory), ylim=c(-400, 400), xlim=c(-10000, 10000),
            zlim=c(-10000, 10000))
})
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()

Looks like this one started out in a strange place.

4.2.2 Muon with many turns

gsNDFSum %>% filter(volCategory == 'Arc', totalDist_meters > 870)

Let’s try event 1509

gsNDF %>% filter(eventEntry == 1509) -> gsNDF1509
gsNDF1509

gsNDF1509 %>% group_by(eventEntry, volCategory) %>% summarize(
                                                totalDist_meters=round(sum(stepLength)/1000, 2), 
                                                totalELoss_MeV=round(sum(totalEnergyDeposit), 1),
                                                nSteps = n()) -> gsNDF1509Sum
gsNDF1509Sum

Let’s plot it!

gsNDF1509 %>% mutate(volCategory = as.factor(as.character(volCategory))) -> gsNDF1509

# lengend3d looks terrible - so do a regular legend
plot(0, type='n', axes=FALSE, ann=FALSE)
with(gsNDF1509,
  legend("bottom", legend=unique(volCategory), col=as.integer(unique(volCategory)), pch=19)
)

clear3d()
with(gsNDF1509,{
     plot3d(x=x, y=y, z=z, type='p', col = as.integer(volCategory), ylim=c(-400, 400), xlim=c(-10000, 10000),
            zlim=c(-10000, 10000))
})
plot3d(ring, add=T, alpha=0.2)
view3d(phi=90, theta=-90)
rglwidget()

Looks like a nice stored muon!! Can we see CBO?

gsNDF1509 %>% mutate(r = sqrt(x*x+z*z)) -> gsNDF1509

p1 <- gsNDF1509 %>% filter(volCategory %in% c('RngTrkPlane')) %>% ggplot(aes(x=globalTime, y=r)) + geom_line() +
  labs(x='Global time (ns)', y='Radius (mm)', title='Radial CBO', subtitle="Measured by Ring Tracking Planes")
#
p2 <- gsNDF1509 %>% filter(volCategory %in% c('Arc','RngTrkPlane')) %>% ggplot(aes(x=globalTime, y=r)) + 
  geom_line() + 
  labs(x='Global time (ns)', y='Radius (mm)', subtitle="Measured by Arc hits and Ring Tracking Planes")
#
grid.arrange(p1, p2)

We can even see the kick!

Here we’ll plot the points too.

p1 <- gsNDF1509 %>% filter(volCategory %in% c('RngTrkPlane')) %>% ggplot(aes(x=globalTime, y=y)) + geom_line() +
  geom_point() +
  labs(x='Global time (ns)', y='y (mm)', title='Vertical CBO', subtitle="Measured by Ring Tracking Planes")
#
p2 <- gsNDF1509 %>% filter(volCategory %in% c('Arc','RngTrkPlane')) %>% ggplot(aes(x=globalTime, y=y)) + 
  geom_line() + geom_point() + 
  labs(x='Global time (ns)', y='y (mm)', subtitle="Measured by Arc hits and Ring Tracking Planes")
#
grid.arrange(p1, p2)